#This program solves the Navier Stokes equation. It uses the method of
#constructed solutions to test for convergence, which should be 2nd order 

#from __future__ import division
import numpy as np
import numpy.fft as fft
from numpy import pi 
from pyfft.cl import Plan
import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.elementwise import ElementwiseKernel
import time
import sys
from pyvisfile.vtk import write_structured_grid
np.set_printoptions(threshold=sys.maxint)
import cPickle
from warnings import filterwarnings
filterwarnings("ignore", module="pyopencl.cache", lineno=336)
filterwarnings("ignore", module="pyopencl", lineno=163)
TIMING = 0 
WGPSIZE_1 = 8 
WGPSIZE_2 = 8 
b_WGPSIZE = 4 
PR_INT = 0
start_size = 5
end_size = 5 

def spread_routine(q, firstn1, nextn1, firstn2, nextn2, F_1, F_2, Mu, X_1, X_2, N, dx, dtheta, f_1, f_2, mu_1, mu_2):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            spread_kernel(q, (N/4, N/4), (WGPSIZE_1, WGPSIZE_2), firstn1.data, nextn1.data, firstn2.data, nextn2.data, F_1.data, F_2.data, Mu.data, X_1.data, X_2.data, N, dx, dtheta, shift0, shift1, f_1.data, f_2.data, mu_1.data, mu_2.data)

def update_routine(q, firstn1, nextn1, firstn2, nextn2, X_1, X_2, N, dx):
    for shift0 in range(0,4):
        for shift1 in range(0,4):
            updaten_kernel(q, (N/4, N/4), (WGPSIZE_1, WGPSIZE_2), firstn1.data, nextn1.data, firstn2.data, nextn2.data, X_1.data, X_2.data, N, dx, shift0, shift1)

def calc_transf(dx):
    N = np.int(np.round(1/dx,0))
    k = np.arange(0, N/2+1)
    k = np.hstack([k, np.arange(-N/2+1, 0)])
    k1 = np.tile(k, (N, 1))
    k2 = k1.T.copy()
    L_tf = -4./dx**2 * (np.sin(pi*k1/N)**2+np.sin(pi*k2/N)**2) 
    L_tf[0,0] = 1. #avoid dividing by zero
    Nyq_filt = np.ones((N,N))
    Nyq_filt[N/2,:]=0
    Nyq_filt[:,N/2]=0 
    p_filt = np.ones((N,N))
    p_filt[0,0] = 0
    return L_tf, Nyq_filt, p_filt

def calc_inds(N): #calculate shifted indices
    a1 = range(0,N)
    i_1 = np.tile(a1, (N, 1))
    i_2 = i_1.T.copy()
    a2 = np.hstack([range(1, N), 0])
    ip_1 = np.tile(a2, (N,1))
    ip_2 = ip_1.T.copy()
    a3 = np.hstack([N-1, range(0, N-1)])
    im_1 = np.tile(a3, (N,1))
    im_2 = im_1.T.copy()
    return i_1, i_2, ip_1, ip_2, im_1, im_2

def solve(t_step):
    global dx, dt, mu_0, rho, b_q, dtheta
    N = np.int(np.round(1/dx,0))
 
    imag_gpu(u_1_gpu, u_1_tf_gpu)
    imag_gpu(u_2_gpu, u_2_tf_gpu)

    plan.execute(u_1_tf_gpu.data)
    plan.execute(u_2_tf_gpu.data)
    
    calc_lagF_gpu(X_1_gpu, X_2_gpu, rr_lag_gpu, ip_s_lag_gpu, im_s_lag_gpu, ds, F_1_gpu, F_2_gpu)

    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)
    zero_gpu(force_1_gpu)
    zero_gpu(force_2_gpu)
    zero_gpu(mu_1_gpu)
    zero_gpu(mu_2_gpu)

    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, Mu_gpu, X_1_gpu, X_2_gpu, N, dx, dtheta, force_1_gpu, force_2_gpu, mu_1_gpu, mu_2_gpu)

    add_gpu(f_1_gpu, force_1_gpu, 1, f_1_gpu) 
    add_gpu(f_2_gpu, force_2_gpu, 1, f_2_gpu) 
    ##########################################################solve for u_tilde
    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, dx, N, S_1_gpu.data, S_2_gpu.data)
    add_gpu(f_1_gpu, S_1_gpu, -rho, f_1_gpu)
    add_gpu(f_2_gpu, S_2_gpu, -rho, f_2_gpu)

    vc_laplace_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_gpu.data, u_2_gpu.data, mu_1_gpu.data, mu_2_gpu.data, dx, N, V_1_gpu.data, V_2_gpu.data) 

    add_gpu(f_1_gpu, V_1_gpu, 1.0, f_1_gpu) 
    add_gpu(f_2_gpu, V_2_gpu, 1.0, f_2_gpu) 
    
    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)

    div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), f_1_gpu.data, f_2_gpu.data, dx, N, w_gpu.data)
    imag_gpu(w_gpu, w_tf_gpu)
    plan.execute(w_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(w_tf_gpu, p_filt_gpu, L_tf_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)

    grad_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), p_gpu.data, dx, N, dp_1_gpu.data, dp_2_gpu.data)
    imag_gpu(dp_1_gpu, dp_1_tf_gpu)
    imag_gpu(dp_2_gpu, dp_2_tf_gpu)
    plan.execute(dp_1_tf_gpu.data)
    plan.execute(dp_2_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tf_gpu.data, u_2_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, dp_1_tf_gpu.data, dp_2_tf_gpu.data, dt, rho, mu_0, N, u_1_tf_tilde_gpu.data, u_2_tf_tilde_gpu.data)

    Nyq_gpu(u_1_tf_tilde_gpu, Nyq_filt_gpu, u_1_tf_tilde_gpu)
    Nyq_gpu(u_2_tf_tilde_gpu, Nyq_filt_gpu, u_2_tf_tilde_gpu)

    plan.execute(u_1_tf_tilde_gpu.data, inverse=True)
    plan.execute(u_2_tf_tilde_gpu.data, inverse=True)
    
    real_gpu(u_1_tf_tilde_gpu, u_1_tilde_gpu)
    real_gpu(u_2_tf_tilde_gpu, u_2_tilde_gpu)
    ############################################################

    ############################################################solve for u_soln
    interp_kernel(queue, (b_q, 1), (b_WGPSIZE, 1), u_1_gpu.data, u_2_gpu.data, X_1_gpu.data, X_2_gpu.data, N, dx, U_1_gpu.data, U_2_gpu.data)
    add_gpu(X_1_gpu, U_1_gpu, dt, X_1_tilde_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt, X_2_tilde_gpu)

    calc_lagF_gpu(X_1_tilde_gpu, X_2_tilde_gpu, rr_lag_gpu, ip_s_lag_gpu, im_s_lag_gpu, ds, F_1_gpu, F_2_gpu)

    zero_gpu(force_1_tilde_gpu)
    zero_gpu(force_2_tilde_gpu)
    zero_gpu(f_1_gpu)
    zero_gpu(f_2_gpu)
    zero_gpu(mu_1_gpu)
    zero_gpu(mu_2_gpu)

    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx)
    spread_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, F_1_gpu, F_2_gpu, Mu_gpu, X_1_tilde_gpu, X_2_tilde_gpu, N, dx, dtheta, force_1_tilde_gpu, force_2_tilde_gpu, mu_1_gpu, mu_2_gpu)

    add_gpu(force_1_gpu, force_1_tilde_gpu, 1, force_1_gpu)
    add_gpu(force_2_gpu, force_2_tilde_gpu, 1, force_2_gpu)

    ss_conv_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, dx, N, S_1_tilde_gpu.data, S_2_tilde_gpu.data)
    add_gpu(S_1_gpu, S_1_tilde_gpu, 1, S_1_gpu)
    add_gpu(S_2_gpu, S_2_tilde_gpu, 1, S_2_gpu)
    add_gpu(f_1_gpu, S_1_gpu, -rho*0.5, f_1_gpu)
    add_gpu(f_2_gpu, S_2_gpu, -rho*0.5, f_2_gpu)

    vc_laplace_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tilde_gpu.data, u_2_tilde_gpu.data, mu_1_gpu.data, mu_2_gpu.data, dx, N, V_1_tilde_gpu.data, V_2_tilde_gpu.data) 

    #print mu_1_gpu.get()
    #print V_1_tilde_gpu.get()
    #print V_2_tilde_gpu.get()
    add_gpu(V_1_gpu, V_1_tilde_gpu, 1, V_1_gpu)
    add_gpu(V_2_gpu, V_2_tilde_gpu, 1, V_2_gpu)
    add_gpu(f_1_gpu, V_1_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, V_2_gpu, 0.5, f_2_gpu)

    add_gpu(f_1_gpu, force_1_gpu, 0.5, f_1_gpu)
    add_gpu(f_2_gpu, force_2_gpu, 0.5, f_2_gpu)

    imag_gpu(f_1_gpu, f_1_tf_gpu)
    imag_gpu(f_2_gpu, f_2_tf_gpu)

    plan.execute(f_1_tf_gpu.data)
    plan.execute(f_2_tf_gpu.data)

    div_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), f_1_gpu.data, f_2_gpu.data, dx, N, w_gpu.data)
    imag_gpu(w_gpu, w_tf_gpu)
    plan.execute(w_tf_gpu.data)

    L_adj_2_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    calc_p_gpu(w_tf_gpu, p_filt_gpu, L_tf_gpu, p_tf_gpu)
    plan.execute(p_tf_gpu.data, inverse=True)
    real_gpu(p_tf_gpu, p_gpu)

    grad_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), p_gpu.data, dx, N, dp_1_gpu.data, dp_2_gpu.data)
    imag_gpu(dp_1_gpu, dp_1_tf_gpu)
    imag_gpu(dp_2_gpu, dp_2_tf_gpu)
    plan.execute(dp_1_tf_gpu.data)
    plan.execute(dp_2_tf_gpu.data)

    L_adj_1_kernel(queue, (1,1), (1,1), L_tf_gpu.data)
    solve_kernel(queue, (N,N), (WGPSIZE_1, WGPSIZE_2), u_1_tf_gpu.data, u_2_tf_gpu.data, L_tf_gpu.data, f_1_tf_gpu.data, f_2_tf_gpu.data, dp_1_tf_gpu.data, dp_2_tf_gpu.data, dt, rho, mu_0, N, u_1_tf_tilde_gpu.data, u_2_tf_tilde_gpu.data)

    Nyq_gpu(u_1_tf_tilde_gpu, Nyq_filt_gpu, u_1_tf_tilde_gpu)
    Nyq_gpu(u_2_tf_tilde_gpu, Nyq_filt_gpu, u_2_tf_tilde_gpu)

    plan.execute(u_1_tf_tilde_gpu.data, inverse=True)
    plan.execute(u_2_tf_tilde_gpu.data, inverse=True)
    
    real_gpu(u_1_tf_tilde_gpu, u_1_gpu)
    real_gpu(u_2_tf_tilde_gpu, u_2_gpu)

    interp_kernel(queue, (b_q, 1), (b_WGPSIZE, 1), u_1_gpu.data, u_2_gpu.data, X_1_tilde_gpu.data, X_2_tilde_gpu.data, N, dx, U_1_soln_gpu.data, U_2_soln_gpu.data)
    add_gpu(U_1_gpu, U_1_soln_gpu, 1, U_1_gpu)
    add_gpu(U_2_gpu, U_2_soln_gpu, 1, U_2_gpu)
    add_gpu(X_1_gpu, U_1_gpu, dt/2, X_1_gpu)
    add_gpu(X_2_gpu, U_2_gpu, dt/2, X_2_gpu)
    update_routine(queue, firstn1_gpu, nextn1_gpu, firstn2_gpu, nextn2_gpu, X_1_gpu, X_2_gpu, N, dx)
    ##############################################################
    return

def err(u_1, u_2, u_1_o, u_2_o, X_1, X_2, X_1_o, X_2_o, p, p_o):
    global dx, dtheta
    NN_o = u_1_o.shape[0]
    errinf1 = 0.
    errinf2 = 0.
    err11 = 0.
    err12 = 0.
    err21 = 0.
    err22 = 0.

    NN_r_o = 3*NN_o/8
    NN_s_o = 75*NN_o/16
    NN_r = 3*NN_o/4
    NN_s = 75*NN_o/8
    X_1_o = X_1_o.reshape(NN_s_o,NN_r_o)
    X_2_o = X_2_o.reshape(NN_s_o,NN_r_o)
    X_1 = X_1.reshape(NN_s,NN_r)
    X_2 = X_2.reshape(NN_s,NN_r)
    errinf1X = 0.
    errinf2X = 0.
    err11X = 0.
    err12X = 0.
    err21X = 0.
    err22X = 0.

    errinfp = 0
    err1p = 0
    err2p = 0

    for i in range(0,NN_o):
        for j in range(0,NN_o):
            if (np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))>errinf1):
                errinf1 = np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))
            if (np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))>errinf2):
                errinf2 = np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))
            err21 += np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))**2
            err22 += np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))**2
            err11 += np.abs(u_1_o[j,i]-0.5*(u_1[2*j,2*i]+u_1[2*j+1,2*i]))
            err12 += np.abs(u_2_o[j,i]-0.5*(u_2[2*j,2*i]+u_2[2*j,2*i+1]))
            if (np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))>errinfp):
                errinfp = np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))
            err2p += np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1]))**2
            err1p += np.abs(p_o[j,i]-0.25*(p[2*j,2*i]+p[2*j+1,2*i]+p[2*j,2*i+1]+p[2*j+1,2*i+1])) 

    err21 = (dx/2)*np.sqrt(err21)
    err22 = (dx/2)*np.sqrt(err22)
    err11 = (dx**2/4)*err11
    err12 = (dx**2/4)*err12
    err2p = (dx/2)*np.sqrt(err2p)
    err1p = (dx**2/4)*err1p

    for i in range(0,NN_r_o):
        for j in range(0,NN_s_o):
            if (np.abs(X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))>errinf1X):
                errinf1X = np.abs(X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1])) 
            if (np.abs(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))>errinf2X):
                errinf2X = np.abs(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1])) 
            err21X += np.abs(X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1]))**2
            err22X += np.abs(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1]))**2
            err11X += np.abs(X_1_o[j,i]-0.25*(X_1[2*j,2*i]+X_1[2*j+1,2*i]+X_1[2*j,2*i+1]+X_1[2*j+1,2*i+1])) 
            err12X += np.abs(X_2_o[j,i]-0.25*(X_2[2*j,2*i]+X_2[2*j+1,2*i]+X_2[2*j,2*i+1]+X_2[2*j+1,2*i+1])) 

    err21X = np.sqrt(dtheta/4)*np.sqrt(err21X)
    err22X = np.sqrt(dtheta/4)*np.sqrt(err22X)
    err11X = (dtheta/4)*err11X
    err12X = (dtheta/4)*err12X

    return errinf1, errinf2, err11, err12, err21, err22, errinf1X, errinf2X, err11X, err12X, err21X, err22X, errinfp, err1p, err2p 

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
Nyq_gpu = ElementwiseKernel(ctx,
	"float2 *u, float *filt, float2 *f",
	"f[i] = filt[i]*u[i]",
	"Nyq_gpu")

calc_lagF_gpu = ElementwiseKernel(ctx,
	"float *X_1, float *X_2, "
	"float *rr_lag, unsigned *ip_s_lag, unsigned *im_s_lag, " 
	"float ds, "
	"float *F_1, float *F_2",
	"F_1[i] = 1./(ds*ds)*(1.+sin(2*pi*rr_lag[i]-pi/2))"
	   "*(X_1[ip_s_lag[i]]-2.*X_1[i]+X_1[im_s_lag[i]]); "
	"F_2[i] = 1./(ds*ds)*(1.+sin(2*pi*rr_lag[i]-pi/2))"
	   "*(X_2[ip_s_lag[i]]-2.*X_2[i]+X_2[im_s_lag[i]])",
        #"F_1[i] = 0.;"
        #"F_2[i] = 0.",
	"calc_lagF_gpu",
	preamble="""
	#define pi (float)3.141592653589793
	""")

complex_mult_gpu = ElementwiseKernel(ctx,
	"float2 *u_1, float2 *u_2, float2 *f",
	"f[i] = complex_mul(u_1[i], u_2[i])",
	"complex_mult_gpu",
	preamble="""
	#define complex_ctr(x, y) (float2)(x, y)
	#define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
	""")

mult_gpu = ElementwiseKernel(ctx,
	"float *u_1, float *u_2, float *f",
	"f[i] = u_1[i]*u_2[i]",
	"mult_gpu")

add_gpu = ElementwiseKernel(ctx,
	"float *u_1, float *u_2, float a, float *f",
	"f[i] = u_1[i]+a*u_2[i]",
	"add_gpu")
    
real_gpu = ElementwiseKernel(ctx,
	"float2 *x, float *z",
	"z[i] = x[i].x",
	"real_gpu")

imag_gpu = ElementwiseKernel(ctx,
	"float *x, float2 *z",
	"z[i].x = x[i];"
	"z[i].y = 0",
	"imag_gpu")

zero_gpu = ElementwiseKernel(ctx,
	"float *x",
	"x[i] = 0",
	"zero_gpu")
 
null_imag_gpu = ElementwiseKernel(ctx,
	"float2 *x, float2 *z",
	"z[i].x = x[i].x;"
	"z[i].y = 0",
	"null_imag_gpu")

calc_p_gpu = ElementwiseKernel(ctx,
	"float2 *w_tf, float *filt, float2 *L_tf, float2 *p_tf",
	"p_tf[i] = filt[i]*complex_div(w_tf[i], L_tf[i])",
	"calc_p_gpu",
	preamble="""
	#define complex_ctr(x, y) (float2)(x, y)
	#define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
	""")

if False:
    solve_gpu = ElementwiseKernel(ctx,
    "float2 *u_1_tf, float2 *u_2_tf, float2 *L_tf,"
    "float2 *f_1_tf, float2 *f_2_tf, float2 *dp_1_tf,"
    "float2 *dp_2_tf, float dt, float rho, float mu,"
    "float2 *u_1_tf_tilde, float2 *u_2_tf_tilde",
    "float2 rhs_1 = complex_mul(1+mu*dt/(2*rho)*L_tf[i], u_1_tf[i])+dt/rho*(f_1_tf[i]-dp_1_tf[i]);"
    "float2 rhs_2 = complex_mul(1+mu*dt/(2*rho)*L_tf[i], u_2_tf[i])+dt/rho*(f_2_tf[i]-dp_2_tf[i]);"
    "float2 lhs = 1-mu*dt/(2*rho)*L_tf[i];"
    "u_1_tf_tilde[i] = complex_div(rhs_1, lhs);"
    "u_2_tf_tilde[i] = complex_div(rhs_2, lhs);",
    "solve_gpu",
    preamble="""
    #define complex_ctr(x, y) (float2)(x, y)
    #define complex_mul(a, b) complex_ctr(mad(-(a).y, (b).y, (a).x * (b).x), mad((a).y, (b).x, (a).x * (b).y))
    #define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
    """)

ib_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    __kernel void interp_save(__global const float *u_1, __global const float *u_2, __global const float *X_1, __global const float *X_2, const unsigned N, const float dx, __global float *U_1, __global float *U_2)
    {
      unsigned k = get_global_id(0);
      float s_1_stag = (X_1[k]+0.5)/dx;
      float s_2_stag = (X_2[k]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float r_2_stag = s_2_stag-i_2_stag;
      float s_1_c = (X_1[k]+0.5-dx/2)/dx;
      float s_2_c = (X_2[k]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      float r_1_c = s_1_c-i_1_c;
      float r_2_c = s_2_c-i_2_c;
      
      int4 ind_1, ind_2;
      ind_1.x = (i_1_stag-1 & N-1);
      ind_1.y = (i_1_stag & N-1);
      ind_1.z = (i_1_stag+1 & N-1);
      ind_1.w = (i_1_stag+2 & N-1);

      ind_2.x = (i_2_c-1 & N-1);
      ind_2.y = (i_2_c & N-1);
      ind_2.z = (i_2_c+1 & N-1);
      ind_2.w = (i_2_c+2 & N-1);

      float16 sb_test;

      sb_test.s0 = u_1[ind_1.x+N*ind_2.x];
      sb_test.s1 = u_1[ind_1.y+N*ind_2.x];
      sb_test.s2 = u_1[ind_1.z+N*ind_2.x];
      sb_test.s3 = u_1[ind_1.w+N*ind_2.x];
      sb_test.s4 = u_1[ind_1.x+N*ind_2.y];
      sb_test.s5 = u_1[ind_1.y+N*ind_2.y];
      sb_test.s6 = u_1[ind_1.z+N*ind_2.y];
      sb_test.s7 = u_1[ind_1.w+N*ind_2.y];
      sb_test.s8 = u_1[ind_1.x+N*ind_2.z];
      sb_test.s9 = u_1[ind_1.y+N*ind_2.z];
      sb_test.sa = u_1[ind_1.z+N*ind_2.z];
      sb_test.sb = u_1[ind_1.w+N*ind_2.z];
      sb_test.sc = u_1[ind_1.x+N*ind_2.w];
      sb_test.sd = u_1[ind_1.y+N*ind_2.w];
      sb_test.se = u_1[ind_1.z+N*ind_2.w];
      sb_test.sf = u_1[ind_1.w+N*ind_2.w];

      float q1 = sqrt(1+4.0*r_1_stag*(1-r_1_stag)); 
      float q2 = sqrt(1+4.0*r_2_c*(1-r_2_c));
      float4 v1, v2;

      v1.x = 3-2*r_1_stag-q1; 
      v1.y = 3-2*r_1_stag+q1; 
      v1.z = 1+2*r_1_stag+q1; 
      v1.w = 1+2*r_1_stag-q1; 
      v1 = v1/8;

      v2.x = 3-2*r_2_c-q2; 
      v2.y = 3-2*r_2_c+q2; 
      v2.z = 1+2*r_2_c+q2; 
      v2.w = 1+2*r_2_c-q2; 
      v2 = v2/8;
      
      float16 w;
      w.s0 = v1.x*v2.x;
      w.s1 = v1.y*v2.x;
      w.s2 = v1.z*v2.x;
      w.s3 = v1.w*v2.x;
      w.s4 = v1.x*v2.y;
      w.s5 = v1.y*v2.y;
      w.s6 = v1.z*v2.y;
      w.s7 = v1.w*v2.y;
      w.s8 = v1.x*v2.z;
      w.s9 = v1.y*v2.z;
      w.sa = v1.z*v2.z;
      w.sb = v1.w*v2.z;
      w.sc = v1.x*v2.w;
      w.sd = v1.y*v2.w;
      w.se = v1.z*v2.w;
      w.sf = v1.w*v2.w;
      
      w = w*sb_test;
      U_1[k] = w.s0+w.s1+w.s2+w.s3
	       +w.s4+w.s5+w.s6+w.s7
	       +w.s8+w.s9+w.sa+w.sb
	       +w.sc+w.sd+w.se+w.sf;

      ind_1.x = (i_1_c-1 & N-1);
      ind_1.y = (i_1_c & N-1);
      ind_1.z = (i_1_c+1 & N-1);
      ind_1.w = (i_1_c+2 & N-1);

      ind_2.x = (i_2_stag-1 & N-1);
      ind_2.y = (i_2_stag & N-1);
      ind_2.z = (i_2_stag+1 & N-1);
      ind_2.w = (i_2_stag+2 & N-1);

      sb_test.s0 = u_2[ind_1.x+N*ind_2.x];
      sb_test.s1 = u_2[ind_1.y+N*ind_2.x];
      sb_test.s2 = u_2[ind_1.z+N*ind_2.x];
      sb_test.s3 = u_2[ind_1.w+N*ind_2.x];
      sb_test.s4 = u_2[ind_1.x+N*ind_2.y];
      sb_test.s5 = u_2[ind_1.y+N*ind_2.y];
      sb_test.s6 = u_2[ind_1.z+N*ind_2.y];
      sb_test.s7 = u_2[ind_1.w+N*ind_2.y];
      sb_test.s8 = u_2[ind_1.x+N*ind_2.z];
      sb_test.s9 = u_2[ind_1.y+N*ind_2.z];
      sb_test.sa = u_2[ind_1.z+N*ind_2.z];
      sb_test.sb = u_2[ind_1.w+N*ind_2.z];
      sb_test.sc = u_2[ind_1.x+N*ind_2.w];
      sb_test.sd = u_2[ind_1.y+N*ind_2.w];
      sb_test.se = u_2[ind_1.z+N*ind_2.w];
      sb_test.sf = u_2[ind_1.w+N*ind_2.w];

      q1 = sqrt(1+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1+4.0*r_2_stag*(1-r_2_stag));

      v1.x = 3-2*r_1_c-q1; 
      v1.y = 3-2*r_1_c+q1; 
      v1.z = 1+2*r_1_c+q1; 
      v1.w = 1+2*r_1_c-q1; 
      v1 = v1/8;

      v2.x = 3-2*r_2_stag-q2; 
      v2.y = 3-2*r_2_stag+q2; 
      v2.z = 1+2*r_2_stag+q2; 
      v2.w = 1+2*r_2_stag-q2; 
      v2 = v2/8;
      
      w.s0 = v1.x*v2.x;
      w.s1 = v1.y*v2.x;
      w.s2 = v1.z*v2.x;
      w.s3 = v1.w*v2.x;
      w.s4 = v1.x*v2.y;
      w.s5 = v1.y*v2.y;
      w.s6 = v1.z*v2.y;
      w.s7 = v1.w*v2.y;
      w.s8 = v1.x*v2.z;
      w.s9 = v1.y*v2.z;
      w.sa = v1.z*v2.z;
      w.sb = v1.w*v2.z;
      w.sc = v1.x*v2.w;
      w.sd = v1.y*v2.w;
      w.se = v1.z*v2.w;
      w.sf = v1.w*v2.w;
      
      w = w*sb_test;
      U_2[k] = w.s0+w.s1+w.s2+w.s3
	       +w.s4+w.s5+w.s6+w.s7
	       +w.s8+w.s9+w.sa+w.sb
	       +w.sc+w.sd+w.se+w.sf;
    }
    __kernel void spread_save(__global const int *firstn, __global const int *nextn, __global const float *F_1, __global const float *F_2, __global const float *X_1, __global const float *X_2, const unsigned N, const float dx, const float dtheta, const unsigned shift0, const unsigned shift1, __global float *f_1, __global float *f_2)
    {
      unsigned i = 8*get_global_id(0)+shift0;
      unsigned j = 8*get_global_id(1)+shift1;
      //unsigned i = 4*get_global_id(0)+shift0;
      //unsigned j = 4*get_global_id(1)+shift1;
      float c = dtheta/(dx*dx);

      int k = firstn[i+N*j]; 
      while(k != -1)
      {
      float s_1_stag = (X_1[k]+0.5)/dx;
      float s_2_stag = (X_2[k]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float r_2_stag = s_2_stag-i_2_stag;
      float s_1_c = (X_1[k]+0.5-dx/2)/dx;
      float s_2_c = (X_2[k]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      float r_1_c = s_1_c-i_1_c;
      float r_2_c = s_2_c-i_2_c;
      
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;
      ind_1x = (i_1_stag-1 & N-1);
      ind_1y = (i_1_stag & N-1);
      ind_1z = (i_1_stag+1 & N-1);
      ind_1w = (i_1_stag+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      float q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      float q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      float v1x, v1y, v1z, v1w;
      float v2x, v2y, v2z, v2w;

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      float ws0, ws1, ws2, ws3;
      float ws4, ws5, ws6, ws7;
      float ws8, ws9, wsa, wsb;
      float wsc, wsd, wse, wsf;
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      float my_F_1 = F_1[k];
      f_1[ind_1x+N*ind_2x] += c*my_F_1*ws0;
      f_1[ind_1y+N*ind_2x] += c*my_F_1*ws1;
      f_1[ind_1z+N*ind_2x] += c*my_F_1*ws2;
      f_1[ind_1w+N*ind_2x] += c*my_F_1*ws3;
      f_1[ind_1x+N*ind_2y] += c*my_F_1*ws4;
      f_1[ind_1y+N*ind_2y] += c*my_F_1*ws5;
      f_1[ind_1z+N*ind_2y] += c*my_F_1*ws6;
      f_1[ind_1w+N*ind_2y] += c*my_F_1*ws7;
      f_1[ind_1x+N*ind_2z] += c*my_F_1*ws8;
      f_1[ind_1y+N*ind_2z] += c*my_F_1*ws9;
      f_1[ind_1z+N*ind_2z] += c*my_F_1*wsa;
      f_1[ind_1w+N*ind_2z] += c*my_F_1*wsb;
      f_1[ind_1x+N*ind_2w] += c*my_F_1*wsc;
      f_1[ind_1y+N*ind_2w] += c*my_F_1*wsd;
      f_1[ind_1z+N*ind_2w] += c*my_F_1*wse;
      f_1[ind_1w+N*ind_2w] += c*my_F_1*wsf;

      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_stag-1 & N-1);
      ind_2y = (i_2_stag & N-1);
      ind_2z = (i_2_stag+1 & N-1);
      ind_2w = (i_2_stag+2 & N-1);

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      float my_F_2 = F_2[k];
      f_2[ind_1x+N*ind_2x] += c*my_F_2*ws0;
      f_2[ind_1y+N*ind_2x] += c*my_F_2*ws1;
      f_2[ind_1z+N*ind_2x] += c*my_F_2*ws2;
      f_2[ind_1w+N*ind_2x] += c*my_F_2*ws3;
      f_2[ind_1x+N*ind_2y] += c*my_F_2*ws4;
      f_2[ind_1y+N*ind_2y] += c*my_F_2*ws5;
      f_2[ind_1z+N*ind_2y] += c*my_F_2*ws6;
      f_2[ind_1w+N*ind_2y] += c*my_F_2*ws7;
      f_2[ind_1x+N*ind_2z] += c*my_F_2*ws8;
      f_2[ind_1y+N*ind_2z] += c*my_F_2*ws9;
      f_2[ind_1z+N*ind_2z] += c*my_F_2*wsa;
      f_2[ind_1w+N*ind_2z] += c*my_F_2*wsb;
      f_2[ind_1x+N*ind_2w] += c*my_F_2*wsc;
      f_2[ind_1y+N*ind_2w] += c*my_F_2*wsd;
      f_2[ind_1z+N*ind_2w] += c*my_F_2*wse;
      f_2[ind_1w+N*ind_2w] += c*my_F_2*wsf;
      k = nextn[k];
      }
    }
    __kernel void interp(__global const float *u_1, __global const float *u_2, __global const float *X_1, __global const float *X_2, const unsigned N, const float dx, __global float *U_1, __global float *U_2)
    {
      unsigned k = get_global_id(0);
      float s_1_stag = (X_1[k]+0.5)/dx;
      float s_2_stag = (X_2[k]+0.5)/dx;       
      int i_1_stag = convert_int(floor(s_1_stag));
      int i_2_stag = convert_int(floor(s_2_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float r_2_stag = s_2_stag-i_2_stag;
      float s_1_c = (X_1[k]+0.5-dx/2)/dx;
      float s_2_c = (X_2[k]+0.5-dx/2)/dx;        
      int i_1_c = convert_int(floor(s_1_c));
      int i_2_c = convert_int(floor(s_2_c));
      float r_1_c = s_1_c-i_1_c;
      float r_2_c = s_2_c-i_2_c;
      
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;
      ind_1x = (i_1_stag-1 & N-1);
      ind_1y = (i_1_stag & N-1);
      ind_1z = (i_1_stag+1 & N-1);
      ind_1w = (i_1_stag+2 & N-1);

      ind_2x = (i_2_c-1 & N-1);
      ind_2y = (i_2_c & N-1);
      ind_2z = (i_2_c+1 & N-1);
      ind_2w = (i_2_c+2 & N-1);

      float sb_tests0, sb_tests1, sb_tests2, sb_tests3;
      float sb_tests4, sb_tests5, sb_tests6, sb_tests7;
      float sb_tests8, sb_tests9, sb_testsa, sb_testsb;
      float sb_testsc, sb_testsd, sb_testse, sb_testsf;

      sb_tests0 = u_1[ind_1x+N*ind_2x];
      sb_tests1 = u_1[ind_1y+N*ind_2x];
      sb_tests2 = u_1[ind_1z+N*ind_2x];
      sb_tests3 = u_1[ind_1w+N*ind_2x];
      sb_tests4 = u_1[ind_1x+N*ind_2y];
      sb_tests5 = u_1[ind_1y+N*ind_2y];
      sb_tests6 = u_1[ind_1z+N*ind_2y];
      sb_tests7 = u_1[ind_1w+N*ind_2y];
      sb_tests8 = u_1[ind_1x+N*ind_2z];
      sb_tests9 = u_1[ind_1y+N*ind_2z];
      sb_testsa = u_1[ind_1z+N*ind_2z];
      sb_testsb = u_1[ind_1w+N*ind_2z];
      sb_testsc = u_1[ind_1x+N*ind_2w];
      sb_testsd = u_1[ind_1y+N*ind_2w];
      sb_testse = u_1[ind_1z+N*ind_2w];
      sb_testsf = u_1[ind_1w+N*ind_2w];

      float q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      float q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));
      float v1x, v1y, v1z, v1w;
      float v2x, v2y, v2z, v2w;

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      float ws0, ws1, ws2, ws3;
      float ws4, ws5, ws6, ws7;
      float ws8, ws9, wsa, wsb;
      float wsc, wsd, wse, wsf;
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      ws0 = ws0*sb_tests0;
      ws1 = ws1*sb_tests1;
      ws2 = ws2*sb_tests2;
      ws3 = ws3*sb_tests3;
      ws4 = ws4*sb_tests4;
      ws5 = ws5*sb_tests5;
      ws6 = ws6*sb_tests6;
      ws7 = ws7*sb_tests7;
      ws8 = ws8*sb_tests8;
      ws9 = ws9*sb_tests9;
      wsa = wsa*sb_testsa;
      wsb = wsb*sb_testsb;
      wsc = wsc*sb_testsc;
      wsd = wsd*sb_testsd;
      wse = wse*sb_testse;
      wsf = wsf*sb_testsf;

      U_1[k] = ws0+ws1+ws2+ws3
	       +ws4+ws5+ws6+ws7
	       +ws8+ws9+wsa+wsb
	       +wsc+wsd+wse+wsf;

      ind_1x = (i_1_c-1 & N-1);
      ind_1y = (i_1_c & N-1);
      ind_1z = (i_1_c+1 & N-1);
      ind_1w = (i_1_c+2 & N-1);

      ind_2x = (i_2_stag-1 & N-1);
      ind_2y = (i_2_stag & N-1);
      ind_2z = (i_2_stag+1 & N-1);
      ind_2w = (i_2_stag+2 & N-1);

      sb_tests0 = u_2[ind_1x+N*ind_2x];
      sb_tests1 = u_2[ind_1y+N*ind_2x];
      sb_tests2 = u_2[ind_1z+N*ind_2x];
      sb_tests3 = u_2[ind_1w+N*ind_2x];
      sb_tests4 = u_2[ind_1x+N*ind_2y];
      sb_tests5 = u_2[ind_1y+N*ind_2y];
      sb_tests6 = u_2[ind_1z+N*ind_2y];
      sb_tests7 = u_2[ind_1w+N*ind_2y];
      sb_tests8 = u_2[ind_1x+N*ind_2z];
      sb_tests9 = u_2[ind_1y+N*ind_2z];
      sb_testsa = u_2[ind_1z+N*ind_2z];
      sb_testsb = u_2[ind_1w+N*ind_2z];
      sb_testsc = u_2[ind_1x+N*ind_2w];
      sb_testsd = u_2[ind_1y+N*ind_2w];
      sb_testse = u_2[ind_1z+N*ind_2w];
      sb_testsf = u_2[ind_1w+N*ind_2w];

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      ws0 = ws0*sb_tests0;
      ws1 = ws1*sb_tests1;
      ws2 = ws2*sb_tests2;
      ws3 = ws3*sb_tests3;
      ws4 = ws4*sb_tests4;
      ws5 = ws5*sb_tests5;
      ws6 = ws6*sb_tests6;
      ws7 = ws7*sb_tests7;
      ws8 = ws8*sb_tests8;
      ws9 = ws9*sb_tests9;
      wsa = wsa*sb_testsa;
      wsb = wsb*sb_testsb;
      wsc = wsc*sb_testsc;
      wsd = wsd*sb_testsd;
      wse = wse*sb_testse;
      wsf = wsf*sb_testsf;

      U_2[k] = ws0+ws1+ws2+ws3
	       +ws4+ws5+ws6+ws7
	       +ws8+ws9+wsa+wsb
	       +wsc+wsd+wse+wsf;
    }
    __kernel void spread(__global const int *firstn1, __global const int *nextn1, __global const int *firstn2, __global const int *nextn2, __global const float *F_1, __global const float *F_2, __global const float *Mu, __global const float *X_1, __global const float *X_2, const unsigned N, const float dx, const float dtheta, const unsigned shift0, const unsigned shift1, __global float *f_1, __global float *f_2, __global float *mu_1, __global float *mu_2)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      float c = dtheta/(dx*dx);
      int ind_1x, ind_1y, ind_1z, ind_1w;
      int ind_2x, ind_2y, ind_2z, ind_2w;

      float f_10, f_11, f_12, f_13;
      float f_14, f_15, f_16, f_17;
      float f_18, f_19, f_1a, f_1b;
      float f_1c, f_1d, f_1e, f_1f;

      float f_20, f_21, f_22, f_23;
      float f_24, f_25, f_26, f_27;
      float f_28, f_29, f_2a, f_2b;
      float f_2c, f_2d, f_2e, f_2f;

      float mu_10, mu_11, mu_12, mu_13;
      float mu_14, mu_15, mu_16, mu_17;
      float mu_18, mu_19, mu_1a, mu_1b;
      float mu_1c, mu_1d, mu_1e, mu_1f;

      float mu_20, mu_21, mu_22, mu_23;
      float mu_24, mu_25, mu_26, mu_27;
      float mu_28, mu_29, mu_2a, mu_2b;
      float mu_2c, mu_2d, mu_2e, mu_2f;

      float v1x, v1y, v1z, v1w;
      float v2x, v2y, v2z, v2w;

      float q1, q2;

      float ws0, ws1, ws2, ws3;
      float ws4, ws5, ws6, ws7;
      float ws8, ws9, wsa, wsb;
      float wsc, wsd, wse, wsf;

      if (k1 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_10=f_1[ind_1x+N*ind_2x];
      f_11=f_1[ind_1y+N*ind_2x];
      f_12=f_1[ind_1z+N*ind_2x];
      f_13=f_1[ind_1w+N*ind_2x];
      f_14=f_1[ind_1x+N*ind_2y];
      f_15=f_1[ind_1y+N*ind_2y];
      f_16=f_1[ind_1z+N*ind_2y];
      f_17=f_1[ind_1w+N*ind_2y];
      f_18=f_1[ind_1x+N*ind_2z];
      f_19=f_1[ind_1y+N*ind_2z];
      f_1a=f_1[ind_1z+N*ind_2z];
      f_1b=f_1[ind_1w+N*ind_2z];
      f_1c=f_1[ind_1x+N*ind_2w];
      f_1d=f_1[ind_1y+N*ind_2w];
      f_1e=f_1[ind_1z+N*ind_2w];
      f_1f=f_1[ind_1w+N*ind_2w];

      mu_10=mu_1[ind_1x+N*ind_2x];
      mu_11=mu_1[ind_1y+N*ind_2x];
      mu_12=mu_1[ind_1z+N*ind_2x];
      mu_13=mu_1[ind_1w+N*ind_2x];
      mu_14=mu_1[ind_1x+N*ind_2y];
      mu_15=mu_1[ind_1y+N*ind_2y];
      mu_16=mu_1[ind_1z+N*ind_2y];
      mu_17=mu_1[ind_1w+N*ind_2y];
      mu_18=mu_1[ind_1x+N*ind_2z];
      mu_19=mu_1[ind_1y+N*ind_2z];
      mu_1a=mu_1[ind_1z+N*ind_2z];
      mu_1b=mu_1[ind_1w+N*ind_2z];
      mu_1c=mu_1[ind_1x+N*ind_2w];
      mu_1d=mu_1[ind_1y+N*ind_2w];
      mu_1e=mu_1[ind_1z+N*ind_2w];
      mu_1f=mu_1[ind_1w+N*ind_2w];
      }

      while(k1 != -1)
      {
      float s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      float r_1_stag = s_1_stag-i_1_stag;
      float s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
      float r_2_c = s_2_c-i_2_c;
      
      q1 = sqrt(1.+4.0*r_1_stag*(1-r_1_stag)); 
      q2 = sqrt(1.+4.0*r_2_c*(1-r_2_c));

      v1x = 3.-2.*r_1_stag-q1; 
      v1y = 3.-2.*r_1_stag+q1; 
      v1z = 1.+2.*r_1_stag+q1; 
      v1w = 1.+2.*r_1_stag-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_c-q2; 
      v2y = 3.-2.*r_2_c+q2; 
      v2z = 1.+2.*r_2_c+q2; 
      v2w = 1.+2.*r_2_c-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      float my_F_1 = F_1[k1];
      f_10 += c*my_F_1*ws0;
      f_11 += c*my_F_1*ws1;
      f_12 += c*my_F_1*ws2;
      f_13 += c*my_F_1*ws3;
      f_14 += c*my_F_1*ws4;
      f_15 += c*my_F_1*ws5;
      f_16 += c*my_F_1*ws6;
      f_17 += c*my_F_1*ws7;
      f_18 += c*my_F_1*ws8;
      f_19 += c*my_F_1*ws9;
      f_1a += c*my_F_1*wsa;
      f_1b += c*my_F_1*wsb;
      f_1c += c*my_F_1*wsc;
      f_1d += c*my_F_1*wsd;
      f_1e += c*my_F_1*wse;
      f_1f += c*my_F_1*wsf;

      float my_Mu = Mu[k1];
      mu_10 += c*my_Mu*ws0;
      mu_11 += c*my_Mu*ws1;
      mu_12 += c*my_Mu*ws2;
      mu_13 += c*my_Mu*ws3;
      mu_14 += c*my_Mu*ws4;
      mu_15 += c*my_Mu*ws5;
      mu_16 += c*my_Mu*ws6;
      mu_17 += c*my_Mu*ws7;
      mu_18 += c*my_Mu*ws8;
      mu_19 += c*my_Mu*ws9;
      mu_1a += c*my_Mu*wsa;
      mu_1b += c*my_Mu*wsb;
      mu_1c += c*my_Mu*wsc;
      mu_1d += c*my_Mu*wsd;
      mu_1e += c*my_Mu*wse;
      mu_1f += c*my_Mu*wsf;

      k1 = nextn1[k1];
      if (k1 == -1){
      f_1[ind_1x+N*ind_2x]=f_10;
      f_1[ind_1y+N*ind_2x]=f_11;
      f_1[ind_1z+N*ind_2x]=f_12;
      f_1[ind_1w+N*ind_2x]=f_13;
      f_1[ind_1x+N*ind_2y]=f_14;
      f_1[ind_1y+N*ind_2y]=f_15;
      f_1[ind_1z+N*ind_2y]=f_16;
      f_1[ind_1w+N*ind_2y]=f_17;
      f_1[ind_1x+N*ind_2z]=f_18;
      f_1[ind_1y+N*ind_2z]=f_19;
      f_1[ind_1z+N*ind_2z]=f_1a;
      f_1[ind_1w+N*ind_2z]=f_1b;
      f_1[ind_1x+N*ind_2w]=f_1c;
      f_1[ind_1y+N*ind_2w]=f_1d;
      f_1[ind_1z+N*ind_2w]=f_1e;
      f_1[ind_1w+N*ind_2w]=f_1f;

      mu_1[ind_1x+N*ind_2x]=mu_10;
      mu_1[ind_1y+N*ind_2x]=mu_11;
      mu_1[ind_1z+N*ind_2x]=mu_12;
      mu_1[ind_1w+N*ind_2x]=mu_13;
      mu_1[ind_1x+N*ind_2y]=mu_14;
      mu_1[ind_1y+N*ind_2y]=mu_15;
      mu_1[ind_1z+N*ind_2y]=mu_16;
      mu_1[ind_1w+N*ind_2y]=mu_17;
      mu_1[ind_1x+N*ind_2z]=mu_18;
      mu_1[ind_1y+N*ind_2z]=mu_19;
      mu_1[ind_1z+N*ind_2z]=mu_1a;
      mu_1[ind_1w+N*ind_2z]=mu_1b;
      mu_1[ind_1x+N*ind_2w]=mu_1c;
      mu_1[ind_1y+N*ind_2w]=mu_1d;
      mu_1[ind_1z+N*ind_2w]=mu_1e;
      mu_1[ind_1w+N*ind_2w]=mu_1f;
      }
      }

      if (k2 !=-1){
      ind_1x = (i-1 & N-1);
      ind_1y = (i & N-1);
      ind_1z = (i+1 & N-1);
      ind_1w = (i+2 & N-1);

      ind_2x = (j-1 & N-1);
      ind_2y = (j & N-1);
      ind_2z = (j+1 & N-1);
      ind_2w = (j+2 & N-1);

      f_20=f_2[ind_1x+N*ind_2x];
      f_21=f_2[ind_1y+N*ind_2x];
      f_22=f_2[ind_1z+N*ind_2x];
      f_23=f_2[ind_1w+N*ind_2x];
      f_24=f_2[ind_1x+N*ind_2y];
      f_25=f_2[ind_1y+N*ind_2y];
      f_26=f_2[ind_1z+N*ind_2y];
      f_27=f_2[ind_1w+N*ind_2y];
      f_28=f_2[ind_1x+N*ind_2z];
      f_29=f_2[ind_1y+N*ind_2z];
      f_2a=f_2[ind_1z+N*ind_2z];
      f_2b=f_2[ind_1w+N*ind_2z];
      f_2c=f_2[ind_1x+N*ind_2w];
      f_2d=f_2[ind_1y+N*ind_2w];
      f_2e=f_2[ind_1z+N*ind_2w];
      f_2f=f_2[ind_1w+N*ind_2w];

      mu_20=mu_2[ind_1x+N*ind_2x];
      mu_21=mu_2[ind_1y+N*ind_2x];
      mu_22=mu_2[ind_1z+N*ind_2x];
      mu_23=mu_2[ind_1w+N*ind_2x];
      mu_24=mu_2[ind_1x+N*ind_2y];
      mu_25=mu_2[ind_1y+N*ind_2y];
      mu_26=mu_2[ind_1z+N*ind_2y];
      mu_27=mu_2[ind_1w+N*ind_2y];
      mu_28=mu_2[ind_1x+N*ind_2z];
      mu_29=mu_2[ind_1y+N*ind_2z];
      mu_2a=mu_2[ind_1z+N*ind_2z];
      mu_2b=mu_2[ind_1w+N*ind_2z];
      mu_2c=mu_2[ind_1x+N*ind_2w];
      mu_2d=mu_2[ind_1y+N*ind_2w];
      mu_2e=mu_2[ind_1z+N*ind_2w];
      mu_2f=mu_2[ind_1w+N*ind_2w];
      }
      
      while (k2 != -1){
      float s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      float r_1_c = s_1_c-i_1_c;
      float s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      float r_2_stag = s_2_stag-i_2_stag;

      q1 = sqrt(1.+4.0*r_1_c*(1-r_1_c)); 
      q2 = sqrt(1.+4.0*r_2_stag*(1-r_2_stag));

      v1x = 3.-2.*r_1_c-q1; 
      v1y = 3.-2.*r_1_c+q1; 
      v1z = 1.+2.*r_1_c+q1; 
      v1w = 1.+2.*r_1_c-q1; 
      v1x = v1x/8.;
      v1y = v1y/8.;
      v1z = v1z/8.;
      v1w = v1w/8.;

      v2x = 3.-2.*r_2_stag-q2; 
      v2y = 3.-2.*r_2_stag+q2; 
      v2z = 1.+2.*r_2_stag+q2; 
      v2w = 1.+2.*r_2_stag-q2; 
      v2x = v2x/8.;
      v2y = v2y/8.;
      v2z = v2z/8.;
      v2w = v2w/8.;
      
      ws0 = v1x*v2x;
      ws1 = v1y*v2x;
      ws2 = v1z*v2x;
      ws3 = v1w*v2x;
      ws4 = v1x*v2y;
      ws5 = v1y*v2y;
      ws6 = v1z*v2y;
      ws7 = v1w*v2y;
      ws8 = v1x*v2z;
      ws9 = v1y*v2z;
      wsa = v1z*v2z;
      wsb = v1w*v2z;
      wsc = v1x*v2w;
      wsd = v1y*v2w;
      wse = v1z*v2w;
      wsf = v1w*v2w;
      
      float my_F_2 = F_2[k2];
      f_20 += c*my_F_2*ws0;
      f_21 += c*my_F_2*ws1;
      f_22 += c*my_F_2*ws2;
      f_23 += c*my_F_2*ws3;
      f_24 += c*my_F_2*ws4;
      f_25 += c*my_F_2*ws5;
      f_26 += c*my_F_2*ws6;
      f_27 += c*my_F_2*ws7;
      f_28 += c*my_F_2*ws8;
      f_29 += c*my_F_2*ws9;
      f_2a += c*my_F_2*wsa;
      f_2b += c*my_F_2*wsb;
      f_2c += c*my_F_2*wsc;
      f_2d += c*my_F_2*wsd;
      f_2e += c*my_F_2*wse;
      f_2f += c*my_F_2*wsf;

      float my_Mu = Mu[k2];
      mu_20 += c*my_Mu*ws0;
      mu_21 += c*my_Mu*ws1;
      mu_22 += c*my_Mu*ws2;
      mu_23 += c*my_Mu*ws3;
      mu_24 += c*my_Mu*ws4;
      mu_25 += c*my_Mu*ws5;
      mu_26 += c*my_Mu*ws6;
      mu_27 += c*my_Mu*ws7;
      mu_28 += c*my_Mu*ws8;
      mu_29 += c*my_Mu*ws9;
      mu_2a += c*my_Mu*wsa;
      mu_2b += c*my_Mu*wsb;
      mu_2c += c*my_Mu*wsc;
      mu_2d += c*my_Mu*wsd;
      mu_2e += c*my_Mu*wse;
      mu_2f += c*my_Mu*wsf;

      k2 = nextn2[k2];
      if (k2 == -1){
      f_2[ind_1x+N*ind_2x]=f_20;
      f_2[ind_1y+N*ind_2x]=f_21;
      f_2[ind_1z+N*ind_2x]=f_22;
      f_2[ind_1w+N*ind_2x]=f_23;
      f_2[ind_1x+N*ind_2y]=f_24;
      f_2[ind_1y+N*ind_2y]=f_25;
      f_2[ind_1z+N*ind_2y]=f_26;
      f_2[ind_1w+N*ind_2y]=f_27;
      f_2[ind_1x+N*ind_2z]=f_28;
      f_2[ind_1y+N*ind_2z]=f_29;
      f_2[ind_1z+N*ind_2z]=f_2a;
      f_2[ind_1w+N*ind_2z]=f_2b;
      f_2[ind_1x+N*ind_2w]=f_2c;
      f_2[ind_1y+N*ind_2w]=f_2d;
      f_2[ind_1z+N*ind_2w]=f_2e;
      f_2[ind_1w+N*ind_2w]=f_2f;

      mu_2[ind_1x+N*ind_2x]=mu_20;
      mu_2[ind_1y+N*ind_2x]=mu_21;
      mu_2[ind_1z+N*ind_2x]=mu_22;
      mu_2[ind_1w+N*ind_2x]=mu_23;
      mu_2[ind_1x+N*ind_2y]=mu_24;
      mu_2[ind_1y+N*ind_2y]=mu_25;
      mu_2[ind_1z+N*ind_2y]=mu_26;
      mu_2[ind_1w+N*ind_2y]=mu_27;
      mu_2[ind_1x+N*ind_2z]=mu_28;
      mu_2[ind_1y+N*ind_2z]=mu_29;
      mu_2[ind_1z+N*ind_2z]=mu_2a;
      mu_2[ind_1w+N*ind_2z]=mu_2b;
      mu_2[ind_1x+N*ind_2w]=mu_2c;
      mu_2[ind_1y+N*ind_2w]=mu_2d;
      mu_2[ind_1z+N*ind_2w]=mu_2e;
      mu_2[ind_1w+N*ind_2w]=mu_2f;
      }
      }
    }
    __kernel void updaten(__global int *firstn1, __global int *nextn1, __global int *firstn2, __global int *nextn2, __global const float *X_1, __global const float *X_2, const unsigned N, const float dx, const unsigned shift0, const unsigned shift1)
    {
      unsigned i = 4*get_global_id(0)+shift0;
      unsigned j = 4*get_global_id(1)+shift1;

      int k1 = firstn1[i+N*j]; 
      int k2 = firstn2[i+N*j]; 
      if (k1 == -1 && k2 == -1){
      return;
      }

      uint ind_1, ind_2;
      int nextnold;

      int nprev1 = -1;
      while(k1 != -1)
      {
      float s_1_stag = (X_1[k1]+0.5)/dx;
      int i_1_stag = convert_int(floor(s_1_stag));
      float s_2_c = (X_2[k1]+0.5-dx/2)/dx;        
      int i_2_c = convert_int(floor(s_2_c));
      
      ind_1 = (i_1_stag & N-1);
      ind_2 = (i_2_c & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev1 = k1;
      k1 = nextn1[k1];
      }
      else {
      nextnold = nextn1[k1];
      nextn1[k1] = firstn1[ind_1+N*ind_2];
      firstn1[ind_1+N*ind_2] = k1;
      k1 = nextnold;
          if (nprev1 == -1) {
          firstn1[i+N*j] = nextnold;
          }
          else {
          nextn1[nprev1] = nextnold;
          }
      }
      }

      int nprev2 = -1;
      while(k2 != -1)
      {
      float s_2_stag = (X_2[k2]+0.5)/dx;        
      int i_2_stag = convert_int(floor(s_2_stag));
      float s_1_c = (X_1[k2]+0.5-dx/2)/dx;
      int i_1_c = convert_int(floor(s_1_c));
      
      ind_1 = (i_1_c & N-1);
      ind_2 = (i_2_stag & N-1);

      if (ind_1 == i && ind_2 == j){
      nprev2 = k2;
      k2 = nextn2[k2];
      }
      else {
      nextnold = nextn2[k2];
      nextn2[k2] = firstn2[ind_1+N*ind_2];
      firstn2[ind_1+N*ind_2] = k2;
      k2 = nextnold;
          if (nprev2 == -1) {
          firstn2[i+N*j] = nextnold;
          }
          else {
          nextn2[nprev2] = nextnold;
          }
      }
      }
    }
    """).build()

solve_prg = cl.Program(ctx, """
    //#pragma OPENCL EXTENSION cl_amd_printf : enable
    #define complex_ctr(x, y) (float2)(x, y)
    #define complex_mul(a, b) complex_ctr(-(a).y*(b).y+(a).x*(b).x, (a).y*(b).x+(a).x*(b).y)
    #define complex_div(a, b) 1/((b).x*(b).x+(b).y*(b).y)*complex_ctr((a).y*(b).y+(a).x*(b).x,(a).y*(b).x-(a).x*(b).y)
    __kernel void solve(__global const float2 *u_1_tf, __global const float2 *u_2_tf, __global const float2 *L_tf,
    __global const float2 *f_1_tf, __global const float2 *f_2_tf, __global const float2 *dp_1_tf,
    __global const float2 *dp_2_tf, const float dt, const float rho, const float mu, const unsigned N,
    __global float2 *u_1_tf_tilde, __global float2 *u_2_tf_tilde)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      unsigned base = i + N*j;
      float2 m;
      m.x = 1.0+mu*dt/(2.0*rho)*L_tf[base].x; 
      m.y = mu*dt/(2.0*rho)*L_tf[base].y; 
      float2 rhs_1 = complex_mul(m, u_1_tf[base])+dt/rho*(f_1_tf[base]-dp_1_tf[base]);
      float2 rhs_2 = complex_mul(m, u_2_tf[base])+dt/rho*(f_2_tf[base]-dp_2_tf[base]);
      float2 lhs;
      lhs.x = 1.0-mu*dt/(2.0*rho)*L_tf[base].x;
      lhs.y = -mu*dt/(2.0*rho)*L_tf[base].y;
      u_1_tf_tilde[base] = complex_div(rhs_1, lhs);
      u_2_tf_tilde[base] = complex_div(rhs_2, lhs);
    }
    __kernel void vc_laplace(__global const float *u_1,
    __global const float *u_2,
    __global const float *mu_1,
    __global const float *mu_2, const float dx, const unsigned N,
    __global float *V_1, __global float *V_2)
    {
    unsigned i = get_global_id(0);
    unsigned j = get_global_id(1);
    unsigned base = i + N*j;
    unsigned base1 = (i+1 & N-1) + N*j;
    unsigned base2 = (i-1 & N-1) + N*j;
    unsigned base3 = i + N*(j+1 & N-1);
    unsigned base4 = i + N*(j-1 & N-1);
    unsigned base5 = (i-1 & N-1) + N*(j+1 & N-1);
    unsigned base6 = (i+1 & N-1) + N*(j-1 & N-1);

    V_1[base] = 1.0/(dx*dx)*(2.0*(0.5*(mu_1[base]+mu_1[base1])*(u_1[base1]-u_1[base])
               -0.5*(mu_1[base]+mu_1[base2])*(u_1[base]-u_1[base2]))
               +0.5*(mu_1[base3]+mu_1[base])*(u_1[base3]-u_1[base])
               -0.5*(mu_1[base]+mu_1[base4])*(u_1[base]-u_1[base4]) 
               +0.5*(mu_1[base3]+mu_1[base])*(u_2[base3]-u_2[base5])
               -0.5*(mu_1[base]+mu_1[base4])*(u_2[base]-u_2[base2]));

    V_2[base] = 1.0/(dx*dx)*(0.5*(mu_2[base]+mu_2[base1])*(u_2[base1]-u_2[base])
               -0.5*(mu_2[base]+mu_2[base2])*(u_2[base]-u_2[base2])
               +0.5*(mu_2[base]+mu_2[base1])*(u_1[base1]-u_1[base6])
               -0.5*(mu_2[base]+mu_2[base2])*(u_1[base]-u_1[base4])
               +2*(0.5*(mu_2[base]+mu_2[base3])*(u_2[base3]-u_2[base])
               -0.5*(mu_2[base]+mu_2[base4])*(u_2[base]-u_2[base4])));

    }   
    """).build()

prg = cl.Program(ctx, """
    __kernel void L_adj_1(__global float2 *L_tf)
    {
      L_tf[0].x = 0.0;
      L_tf[0].y = 0.0;
    }
    __kernel void L_adj_2(__global float2 *L_tf)
    {
      L_tf[0].x = 1.0;
      L_tf[0].y = 0.0;
    }
    __kernel void div(__global const float *u_1,
    __global const float *u_2,
    const float dx, const unsigned N, __global float *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //f[i,j] = 1/dx*(u_1[i_2[i,j],ip_1[i,j]]-u_1[i,j]+u_2[ip_2[i,j],i_1[i,j]]-u_2[i,j];
      unsigned base = i + N*j;
      unsigned base1 = (i+1 & N-1)+N*j;
      unsigned base2 = i+N*(j+1 & N-1);
      f[base] = 1/dx*(u_1[base1]-u_1[base]+u_2[base2]-u_2[base]);
    }
    __kernel void grad(__global const float *phi,
    const float dx, const unsigned N,
    __global float *g_1, __global float *g_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //grad_1 = 1/dx*(phi[i_2,i_1]-phi[i_2,im_1])
      //grad_2 = 1/dx*(phi[i_2,i_1]-phi[im_2,i_1])
      unsigned base = i+N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = i+N*(j-1 & N-1);
      g_1[base] = 1/dx*(phi[base]-phi[base1]);
      g_2[base] = 1/dx*(phi[base]-phi[base2]);
    }
    __kernel void ss_conv(__global const float *u_1,
    __global const float *u_2, const float dx, const unsigned N,
    __global float *S_1, __global float *S_2)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //Compute S = 0.5*(div(UU) + (U*grad)U).
      unsigned base = i+N*j;
      unsigned base1 = (i-1 & N-1)+N*j;
      unsigned base2 = (i+1 & N-1)+N*j;
      unsigned base3 = i+N*(j-1 & N-1);
      unsigned base4 = i+N*(j+1 & N-1);
      unsigned base5 = (i-1 & N-1)+N*(j+1 & N-1);
      unsigned base6 = (i+1 & N-1)+N*(j-1 & N-1);
      float D_x, D_y, A_x, A_y, D, A;
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]+u_1[base]) -
          (u_1[base]+u_1[base1])*(u_1[base]+u_1[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base5])*(u_1[base4]+u_1[base]) -
          (u_2[base]+u_2[base1])*(u_1[base]+u_1[base3]));
      D = D_x + D_y;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base])*(u_1[base2]-u_1[base]) +
          (u_1[base]+u_1[base1])*(u_1[base]-u_1[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base5])*(u_1[base4]-u_1[base]) +
            (u_2[base]+u_2[base1])*(u_1[base]-u_1[base3]));
      A = A_x + A_y;
      S_1[base] = 0.5*(D+A);
      D_x = (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]+u_2[base]) -
          (u_1[base]+u_1[base3])*(u_2[base]+u_2[base1]));
      D_y = (0.25/dx)*(
          (u_2[base4]+u_2[base])*(u_2[base4]+u_2[base]) -
          (u_2[base]+u_2[base3])*(u_2[base]+u_2[base3]));
      D = D_x + D_y;
      A_x =  (0.25/dx)*(
          (u_1[base2]+u_1[base6])*(u_2[base2]-u_2[base]) +
          (u_1[base]+u_1[base3])*(u_2[base]-u_2[base1]));
      A_y = (0.25/dx)*(
            (u_2[base4]+u_2[base])*(u_2[base4]-u_2[base]) +
            (u_2[base]+u_2[base3])*(u_2[base]-u_2[base3]));
      A = A_x + A_y;
      S_2[base] = 0.5*(D+A);
    }
    __kernel void interpctolr(__global const float *phi, __global const unsigned *i_2,  __global const unsigned *im_1, const unsigned N, __global float *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[i_2,im_1])

      unsigned base = i+N*j;
      unsigned base1 = im_1[base]+N*i_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interpctotb(__global const float *phi, __global const unsigned *i_1, __global const unsigned *im_2, const unsigned N, __global float *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[im_2,i_1])

      unsigned base = i+N*j;
      unsigned base1 = i_1[base]+N*im_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interplrtoc(__global const float *phi, __global const unsigned *i_2, __global const unsigned *ip_1, const unsigned N, __global float *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[i_2,ip_1])

      unsigned base = i+N*j;
      unsigned base1 = ip_1[base]+N*i_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    __kernel void interptbtoc(__global const float *phi, __global const unsigned *i_1,__global const unsigned *ip_2, const unsigned N, __global float *f)
    {
      unsigned i = get_global_id(0);
      unsigned j = get_global_id(1);
      //phi_itp = 0.5*(phi+phi[ip_2,i_1])

      unsigned base = i+N*j;
      unsigned base1 = i_1[base]+N*ip_2[base];
      f[base] = 0.5*(phi[base]+phi[base1]);
    }
    """).build()

L_adj_1_kernel = prg.L_adj_1
L_adj_1_kernel.set_scalar_arg_dtypes([None])

L_adj_2_kernel = prg.L_adj_2
L_adj_2_kernel.set_scalar_arg_dtypes([None])

div_kernel = prg.div
div_kernel.set_scalar_arg_dtypes([None]*2 + [np.float32, np.uint32] + [None])

solve_kernel = solve_prg.solve
solve_kernel.set_scalar_arg_dtypes([None]*7 + [np.float32]*3 + [np.uint32] + [None]*2)

vc_laplace_kernel = solve_prg.vc_laplace
vc_laplace_kernel.set_scalar_arg_dtypes([None]*4 + [np.float32] + [np.uint32] + [None]*2)

grad_kernel = prg.grad
grad_kernel.set_scalar_arg_dtypes([None]*1 + [np.float32, np.uint32] + [None]*2)

ss_conv_kernel = prg.ss_conv
ss_conv_kernel.set_scalar_arg_dtypes([None]*2 + [np.float32, np.uint32] + [None]*2)

interp_kernel = ib_prg.interp
interp_kernel.set_scalar_arg_dtypes([None]*4 + [np.uint32, np.float32] + [None]*2)
spread_kernel = ib_prg.spread
spread_kernel.set_scalar_arg_dtypes([None]*9 + [np.uint32] + [np.float32]*2 + [np.uint32]*2 + [None]*4)
updaten_kernel = ib_prg.updaten
updaten_kernel.set_scalar_arg_dtypes([None]*6 + [np.uint32] + [np.float32] + [np.uint32]*2)

interpctolr_kernel = prg.interpctolr
interpctolr_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interpctotb_kernel = prg.interpctotb
interpctotb_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interplrtoc_kernel = prg.interplrtoc
interplrtoc_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])
interptbtoc_kernel = prg.interptbtoc
interptbtoc_kernel.set_scalar_arg_dtypes([None]*3 + [np.uint32] + [None])

for ii in range(start_size,end_size+1): 
    #set up grid
    N = 2**ii
    dx = 1.0/N
    rho = 1.0

    x_1_stag = np.arange(-0.5, 0.5, dx)
    x_2_stag = np.arange(-0.5, 0.5, dx)
    xx_1_stag, xx_2_stag = np.meshgrid(x_1_stag, x_2_stag)

    x_1_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    x_2_c = np.arange(-0.5+dx/2, 0.5+dx/2, dx)
    xx_1_c, xx_2_c = np.meshgrid(x_1_c, x_2_c)
    
    if PR_INT:
        mesh1 = np.vstack((xx_1_stag, xx_2_c))
        mesh2 = np.vstack((xx_1_c, xx_2_stag))
        mesh3 = np.vstack((xx_1_c, xx_2_c))

    mu_0 = 0.05
    dt = 0.08*dx
    clock_max = np.int(np.round(0.4/dt,0))
    #clock_max = 1
    print clock_max

    #set up Lagrangian parameters
    #N_r = N/2
    #N_s = N/2
    N_r = 3*N/8
    N_s = 75*N/16
    b_q = N_r*N_s
    dr = 1./N_r
    ds = 1./N_s
    dtheta = dr*ds
    a_lag = 0.2
    b_lag = 0.25
    g_lag = 0.0625

    #set up Lagrangian grid
    r_lag = np.arange(dr/2, 1, dr)
    s_lag = np.arange(ds/2, 1, ds)
    rr_lag, ss_lag = np.meshgrid(r_lag, s_lag)
    rr_lag = rr_lag.reshape(b_q,1)
    ss_lag = ss_lag.reshape(b_q,1)

    i_lag = np.arange(0, b_q)
    ip_s_lag = np.remainder(i_lag+N_r, b_q)
    im_s_lag = np.remainder(i_lag-N_r, b_q)

    X_1 = (a_lag+g_lag*(rr_lag-0.5))*np.cos(2*pi*ss_lag)
    X_2 = (b_lag+g_lag*(rr_lag-0.5))*np.sin(2*pi*ss_lag)
    Mu = 0.04*mu_0*(1.+np.sin(2*pi*rr_lag-pi/2))
    #Mu = 0*mu_0*(1.+np.sin(2*pi*rr_lag-pi/2))

    firstn1 = -1*np.ones((N,N))
    nextn1 = -1*np.ones((b_q,1))
    firstn2 = -1*np.ones((N,N))
    nextn2 = -1*np.ones((b_q,1))

    for k in range(0,b_q):
        s_1_stag = (X_1[k]+0.5)/dx
        s_2_stag = (X_2[k]+0.5)/dx        
        i_1_stag = (np.floor(s_1_stag)).astype(np.uint32)
        i_2_stag = (np.floor(s_2_stag)).astype(np.uint32)
        s_1_c = (X_1[k]+0.5-dx/2)/dx;
        s_2_c = (X_2[k]+0.5-dx/2)/dx;        
        i_1_c = (np.remainder(np.floor(s_1_c),N)).astype(np.uint32);
        i_2_c = (np.remainder(np.floor(s_2_c),N)).astype(np.uint32);
        if (firstn1[i_2_c,i_1_stag]==-1):
            firstn1[i_2_c,i_1_stag]=k 
        else:
            nextn1[k]=firstn1[i_2_c,i_1_stag]
            firstn1[i_2_c,i_1_stag]=k
        if (firstn2[i_2_stag,i_1_c]==-1):
            firstn2[i_2_stag,i_1_c]=k 
        else:
            nextn2[k]=firstn2[i_2_stag,i_1_c]
            firstn2[i_2_stag,i_1_c]=k

    L_tf, Nyq_filt, p_filt = calc_transf(dx)
    i_1, i_2, ip_1, ip_2, im_1, im_2 = calc_inds(N)

    #u_1_soln = 2*np.cos(2*pi*xx_1_stag)*np.sin(4*pi*xx_2_c)
    #u_2_soln = -np.sin(2*pi*xx_1_c)*np.cos(4*pi*xx_2_stag)
    u_1_soln = np.zeros((N,N)) #initial conditions
    u_2_soln = np.zeros((N,N))

    #move data to gpu
    plan = Plan((N, N), queue=queue)

    u_1_gpu = cl_array.to_device(ctx, queue, u_1_soln.astype(np.float32))
    u_2_gpu = cl_array.to_device(ctx, queue, u_2_soln.astype(np.float32))
    u_1_tf_gpu = cl_array.to_device(ctx, queue, u_1_soln.astype(np.complex64))
    u_2_tf_gpu = cl_array.to_device(ctx, queue, u_2_soln.astype(np.complex64))
    L_tf_gpu = cl_array.to_device(ctx, queue, L_tf.astype(np.complex64))
    Nyq_filt_gpu = cl_array.to_device(ctx, queue, Nyq_filt.astype(np.float32))
    p_filt_gpu = cl_array.to_device(ctx, queue, p_filt.astype(np.float32))
    firstn1_gpu = cl_array.to_device(ctx, queue, firstn1.astype(np.int32))
    nextn1_gpu = cl_array.to_device(ctx, queue, nextn1.astype(np.int32))
    firstn2_gpu = cl_array.to_device(ctx, queue, firstn2.astype(np.int32))
    nextn2_gpu = cl_array.to_device(ctx, queue, nextn2.astype(np.int32))
    f_1_gpu = cl_array.empty_like(u_1_gpu) 
    f_2_gpu = cl_array.empty_like(u_1_gpu) 
    f_1_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    f_2_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    force_1_gpu = cl_array.empty_like(u_1_gpu) 
    force_2_gpu = cl_array.empty_like(u_1_gpu) 
    mu_1_gpu = cl_array.empty_like(u_1_gpu) 
    mu_2_gpu = cl_array.empty_like(u_1_gpu) 
    force_1_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    force_2_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    p_gpu = cl_array.empty_like(u_1_gpu) 
    p_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    S_1_gpu = cl_array.empty_like(u_1_gpu) 
    S_2_gpu = cl_array.empty_like(u_1_gpu) 
    S_1_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    S_2_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    V_1_gpu = cl_array.empty_like(u_1_gpu) 
    V_2_gpu = cl_array.empty_like(u_1_gpu) 
    V_1_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    V_2_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    w_gpu = cl_array.empty_like(u_1_gpu) 
    w_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    dp_1_gpu = cl_array.empty_like(u_1_gpu) 
    dp_2_gpu = cl_array.empty_like(u_1_gpu) 
    dp_1_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    dp_2_tf_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_1_tf_tilde_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_2_tf_tilde_gpu = cl_array.empty_like(u_1_tf_gpu) 
    u_1_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    u_2_tilde_gpu = cl_array.empty_like(u_1_gpu) 
    X_1_gpu = cl_array.to_device(ctx, queue, X_1.astype(np.float32))
    X_2_gpu = cl_array.to_device(ctx, queue, X_2.astype(np.float32))
    X_1_tilde_gpu = cl_array.to_device(ctx, queue, X_1.astype(np.float32))
    X_2_tilde_gpu = cl_array.to_device(ctx, queue, X_2.astype(np.float32))
    rr_lag_gpu = cl_array.to_device(ctx, queue, rr_lag.astype(np.float32))
    ip_s_lag_gpu = cl_array.to_device(ctx, queue, ip_s_lag.astype(np.uint32))
    im_s_lag_gpu = cl_array.to_device(ctx, queue, im_s_lag.astype(np.uint32))
    F_1_gpu = cl_array.empty_like(X_1_gpu) 
    F_2_gpu = cl_array.empty_like(X_1_gpu) 
    Mu_gpu = cl_array.to_device(ctx, queue, Mu.astype(np.float32)) 
    U_1_gpu = cl_array.empty_like(X_1_gpu) 
    U_2_gpu = cl_array.empty_like(X_1_gpu) 
    U_1_soln_gpu = cl_array.empty_like(X_1_gpu) 
    U_2_soln_gpu = cl_array.empty_like(X_1_gpu) 


    #solve
    queue.finish()
    start_time = time.time()
    for t_step in range(1, clock_max+1):
	solve(t_step)
	if t_step % 50 == 0:
	    print t_step
            #print mu_1_gpu.get().max()
            #print u_1_gpu.get().max()
        if PR_INT:
	    if t_step % PR_INT == 0:
                write_structured_grid("p.vts", mesh3, 
                point_data=[("p", p_gpu.get())])
    queue.finish()
    print 'time spent is', time.time()-start_time, 'seconds'

    u_1_soln = u_1_gpu.get()
    u_2_soln = u_2_gpu.get()
    X_1_soln = X_1_gpu.get()
    X_2_soln = X_2_gpu.get()
    p_soln = p_gpu.get() 
    #print u_1_soln
    #print u_2_soln
    #calculate error
    if ii != start_size: 
        errinf1, errinf2, err11, err12, err21, err22, errinf1X, errinf2X, err11X, err12X, err21X, err22X, errinfp, err1p, err2p = err(u_1_soln, u_2_soln, u_1_o, u_2_o, X_1_soln, X_2_soln, X_1_o, X_2_o, p_soln, p_o)
        #errinf1, errinf2, err21, err22 = err(u_1_soln, u_2_soln, u_1_o, u_2_o)
        print 'u Errinf1 is', errinf1, '\n u Errinf2 is', errinf2
        print 'u Err11 is', err11, '\n u Err12 is', err12
        print 'u Err21 is', err21, '\n u Err22 is', err22
        print 'X Errinf1 is', errinf1X, '\n X Errinf2 is', errinf2X
        print 'X Err11 is', err11X, '\n X Err12 is', err12X
        print 'X Err21 is', err21X, '\n X Err22 is', err22X
        print 'p Errinf is', errinfp
        print 'p Err1 is', err1p
        print 'p Err2 is', err2p
        if ii > start_size+1:
            print 'rinf for u1', -np.log(errinf1/errinf1_o)/np.log(2)
            print 'rinf for u2', -np.log(errinf2/errinf2_o)/np.log(2)
            print 'r1 for u1', -np.log(err11/err11_o)/np.log(2)
            print 'r1 for u2', -np.log(err12/err12_o)/np.log(2)
            print 'r2 for u1', -np.log(err21/err21_o)/np.log(2)
            print 'r2 for u2', -np.log(err22/err22_o)/np.log(2)
            print 'rinf for X1', -np.log(errinf1X/errinf1X_o)/np.log(2)
            print 'rinf for X2', -np.log(errinf2X/errinf2X_o)/np.log(2)
            print 'r1 for X1', -np.log(err11X/err11X_o)/np.log(2)
            print 'r1 for X2', -np.log(err12X/err12X_o)/np.log(2)
            print 'r2 for X1', -np.log(err21X/err21X_o)/np.log(2)
            print 'r2 for X2', -np.log(err22X/err22X_o)/np.log(2)
            print 'rinf for p', -np.log(errinfp/errinfp_o)/np.log(2)
            print 'r1 for p', -np.log(err1p/err1p_o)/np.log(2)
            print 'r2 for p', -np.log(err2p/err2p_o)/np.log(2)
        errinf1_o = errinf1
        errinf2_o = errinf2
        err11_o = err11
        err12_o = err12
        err21_o = err21
        err22_o = err22
        errinf1X_o = errinf1X
        errinf2X_o = errinf2X
        err11X_o = err11X
        err12X_o = err12X
        err21X_o = err21X
        err22X_o = err22X
        errinfp_o = errinfp
        err1p_o = err1p
        err2p_o = err2p
    u_1_o = u_1_soln
    u_2_o = u_2_soln
    X_1_o = X_1_soln
    X_2_o = X_2_soln
    p_o = p_soln

    #rm -f single.pkl
    output = open('single.pkl', 'wb')
    pklist = [u_1_soln, u_2_soln, X_1_soln, X_2_soln, p_soln] 
    cPickle.dump(pklist, output, 2)
    output.close()
